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Abstract 

We discuss the equilibrium of a single collective variable characterizing a finite set of coupled, 
noisy, bistable systems as the noise strength, the size and the coupling parameter are varied. We 
identify distinct regions in parameter space. The results obtained in prior works in the asymptotic 
infinite size limit are significantly different from the finite size results. A procedure to construct 
approximate 1-dimensional Langevin equation is adopted. This equation provides a useful tool to 
understand the collective behavior even in the presence of an external driving force. 



Stochastic resonance (SR) is a phenomenon where the response of a nonlinear dynamical 
system to an external driving is enhanced by the action of noise [l]. Although SR has been 
mainly discussed for simple systems, its analysis has also been extended to complex systems 
with many interacting units 2], [], Q]. Our present study is prompted by recent studies on 



SR in finite arrays J5| 



flfl 



8|. Typically, the noise strength is the parameter varied to 



observe SR effeets In arrays, the coupling strength |3, B, as we., as the system size 
[4], 15| have also been considered as parameters leading to SR. The complexity of the SR 
analysis in arrays can be facilitated by an adequate understanding of the equilibrium state. 
In this work, we carry out a detailed numerical analysis of the equilibrium distribution of the 
collective variable. The reduction of the multidimensional Langevin dynamics to an effective 
1-dimensional Langevin equation greatly simplifies the analysis of the system response to 
weak driving forces 0, Q]. Here we will assess the possibility of such a reduced Langevin 
description. 

Our model consists of a set of N identical bistable units, each of them characterized by 
a variable Xi(t) (i — h. . .,N) satisfying a stochastic evolution equation (in dimensionless 
form) of the type fl, E, Q 

9 N 



x i = x i -x? + —^{x j -x i )+V2DZi{t), (1) 
iV 3=1 

where 9 is a coupling parameter and the term £j (t) represents a white noise with zero average 
and = ~ s )- The set Xi(t) is an iV-dimensional Markovian process. 

We are interested in the properties of a collective variable, S(t) = jjJ2j x j(t), charac- 
terizing the chain as a whole. Even though the set Xi(t) is an N-dimensional Markovian 
process, S(t) in general is not. Consequently, there is no reason why S(t) should sat- 
isfy a Langevin equation. The equilibrium probability density of the collective variable is 
P eq (s) = (5(s — S(t))) eq , where the average is taken with the iV-dimensional equilibrium 
distribution for the Xi(t) process. An exact analytical evaluation of the multidimensional 
integral is, in general, impossible. 

In the limit N — > 00, Desai and Zwanzig [9( carried out an asymptotic analysis of the 
equilibrium density P^(s) based on a saddle point approximation. Their analysis shows 
that P^(s) remains non-Gaussian even if the limit 9 — > is taken at the end. A convenient 
parameter is z = l^gf- The asymptotic P^(s) is not necessarily unique. Indeed, there is 
a critical line in a z vs. 9 diagram such that for z values below the critical line, there is 



a single stable P^(s) which is symmetrical around s = and either bimodal (for 9 < 1) 
or monomodal (for 9 > 1). On the other hand, for any z value above the critical line, 
there are two stable coexisting monomodal distributions centered around values ±s , plus 
one unstable distribution centered at s = 0. Desai and Zwanzig also discuss a Gaussian 
approximation to the general saddle point expression leading to a critical line 9 = 3D. 
Points on this critical line in the z vs. 9 plane are depicted in Fig. (JT^,) and in the D vs. 9 
plane in Fig. (PJ>) (triangles). 

For finite systems, we rely on numerical simulations to obtain information on P eq (s). 
The set of equations in Eq. (CQ) are numerically integrated for very many realizations of the 
noise terms. After a transient period we start gathering data and average over the noise 
realizations to construct histograms estimating P eq (s). We also evaluate the equilibrium 
time correlation function of the collective variable and check that the system has indeed 
relaxed to an equilibrium situation. The H-theorem for finite systems guarantees 11[ that 
the equilibrium distribution function of the iV-dimensional Xi(t) process is unique regardless 
of the initial condition. Thus, it follows that P eq (s) also exists and it is unique. Nonetheless, 
its functional form might depend on the parameter values considered. 

Our numerical findings indicate that there exists a line separating different regions in 
parameter space. In Fig. flT^,) we depict the line for N = 10 (circles), and for N = 50 
(squares) in a 2 vs. 9 plot. For points above the line, P eq (s) has a local minimum at s = 
separating two maxima symmetrically located around zero. Below the line, P eq (s) has a 
single maximum at s = 0. In the D vs. 9 plot in Fig. ([lb), the equilibrium probability density 
for parameter values above the depicted line is always monomodal, with a maximum at s — 0, 
while it is multimodal for points below the line. An example of this exchange of shape as D 
is varied with 9 kept constant can be seen in Fig. ([2]) for 9 — 0.5 and N = 10. The transitions 
among barriers in the iV-dimensional energy surface are induced by the noise. Then, for 
large noise strengths, the random trajectories of Xi(t) have good chances of exploring all 
the attractors quite frequently with numerous jumps over the barriers, leading to a single 
maximum distribution for the global variable. For low noise strengths the bimodality of 
the distribution is to be expected with more scarce jumps over high energy regions. It 
should be pointed out that the shape of P eq (s) does not imply the shape for the one- 
particle equilibrium distribution function /i q (x) obtained by integrating the joint equilibrium 
probability distribution over all the degrees of freedom except one. An example is shown 
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FIG. 1: (Color online)Numerically determined lines separating different regions in parameter space. 
In the z vs. 9 plane (a) P eq (s) is bimodal above the line and monomodal below it. In the z vs. D 
plane (b), P eq (s) is monomodal above the line or multimodal below it. The circles correspond to a 
system with N = 10, while the squares are for N = 50. The coefficient b of the effective Langevin 
equation, Eq. ([2]) does not exist above the line with diamonds. The line with triangles correspond 

n 

to the Gaussian approximation D = 36 in [9] for an infinite size system. 

in Fig. ([3]), where we note that for the parameter values iV = 50, D = 0.15, and 9 = 0.2, 
P eq (s) is monomodal, while /i q (x) is bimodal. 

The features just described for finite systems are in sharp contrast with the behavior found 
in previous works within the N — > oo limit [h]]. The existence of two possible equilibrium 
distributions above the critical line is impossible in finite systems. Also, in the asymptotic 
infinite N limit, low values of z favor a single bimodal distribution for 9 < 1, or monomodal 
for 9 > 1, while for finite values of iV the distribution P eq (s) is monomodal for all 9. On 
the other hand, the two possible single maximum stable equilibrium distributions found for 
large values of z in the asymptotic limit are replaced by a single bimodal distribution in the 
finite size case. 

It should also be pointed out that the lines separating the different regions in either the z 
or D vs. 9 plots are quite insensitive to the N value used, as long as it remains finite. Thus, 
the change in the shape of the global equilibrium distribution does not seem to depend much 
on the system size, as long as it remains finite. 

In the asymptotic limit (N — > oo), Desai and Zwanzig also showed that the dynamics 
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FIG. 2: P eq (s) for a N = 10 units system with 9 = 0.5 and several values of D, obtained from the 
simulations of the full set of equations, Eq. ([!]) 
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FIG. 3: Plots of the equilibrium distribution for the collective variable P eq (s) and the equilibrium 
distribution for one of the subsystems, /^ q (x), for a system with N = 50 units coupled with a 
coupling strength 9 = 0.2 and with D = 0.15. 

could be casted in terms of a truly nonlinear Fokker-Planck equation (i. e., nonlinear in 
the probability distribution) consistent with the bifurcation of the equilibrium probability 
distribution. Further details about the description of the system in the N oo limit 
were discussed later by analytical studies or numerical simulations 10J. There have been 
several attempts to construct effective 1-dimensional Langevin dynamics for the collective 
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variable for finite systems 0, Q]. Although such an equation does not necessarily exist as 
S(t) is not necessarily a Markovian process, when it does, it might be a useful and reliable 
approximation. In 5], Pikovsky et al. used the Gaussian truncation of an infinite hierarchy 
of equations for the cumulants and a slaving principle to construct an effective 1-dimensional 
Langevin equation. Its explicit form is 



2D 

S = aS-bS 3 +J— X (t), (2) 



with (x(t)x( s )) — d{t — s ) an d the coefficients a and b given by 



a = 1 + 0.5(9 - 1) - 0.5^(0 - l) 2 + 12D; 

b = / a (3) 

2 - 9 + J (2 + 9) 2 - 24D 



Unfortunately, the drift term appearing in such equation is a complex number for D values 
below the line labelled as "effective 1-D" in the z-9 plane in Fig. (JTk) (or above the corre- 
sponding line in the D-9 plane in Fig. ((Ib)) and so the effective Langevin equation, Eq. (T5]) 
does not exist for all values in parameter space. 



The method described in [13|] opens up another possibility to construct an effective 1- 
dimensional Langevin equation by numerically estimating the Kramers-Moyal coefficients 
D (n \x,t) for a given data set of a general Markovian process X(t). These coefficients are 
defined as 

DW(x, t) = - lim - ({X(t + r) - x) n > \ x(t)=x . (4) 

As described in [13[| these conditional moments can be estimated for the smallest available 
r by computing histograms. Given that the collective variable S(t) is a Markovian process 
on this scale, i.e., the Markov length is smaller than r, this procedure yields a reliable 
1-dimensional Langevin description of the form 



S = DW( 8 ,t) + y/2DW(s,t) X (t). (5) 

We find that the fourth coefficient is about four orders of magnitude smaller than D^ 2 , 
i. e., it practically vanishes. Hence, the Pawula theorem guarantees that the third and all 
higher coefficients do also vanish. 

In Fig. H we depict our numerically estimated D^(s) for the system in Eq. ([1]) with 
iV = 10, 9 = 0.5 and varying noise strengths. The noise value D c w 0.26 corresponds to 



-1.5 -1 -0.5 0.5 1 1.5 
S 

FIG. 4: (Color online) The drift coefficient D^\s) for the effective Langevin, Eq. ©, for a system 
described by Eq. (pQ) with N = 10, 9 = 0.5 and several D values 

the transition point in the D vs. 9 diagram. The drift coefficients can be fitted with fifth 
degree odd polynomials, in contrast to the result derived by Pikovsky et al. The diffusion 
coefficient is practically constant and satisfies 

= -. (6) 

This relation is also valid for all other parameters we have tested. 

We have also generated histograms for the equilibrium distribution using the effective 
Langevin equation Eq. (jSJ) for parameter values above and below the transition line in Fig. 
([T]b). The results are shown in Fig. ((Sj). Comparing with Fig. (j2j) we see that the effective 
Langevin equation in Eq. (j3J) reproduces quite correctly the results obtained from the full 
set of equations for points above and below the transition line. 

In previous work [7j], we have analyzed the phenomenon of SR for the collective variable of 
an array driven by weak periodic forces. It was demonstrated that when the noise strength 
is varied while keeping 9 and iV constant, the signal-to-noise ratio (SNR) of S(t) shows 
a nonmonotonic behavior with D. Even for weak driving amplitudes, the SNR reaches 
values much larger than those typically observed in single unit systems for the same driving 
forces. This enhancement is largely due to the strong reduction of the fluctuations in the 
driven system with respect to those present in the absence of driving. The multidimensional 
character of the full potential relief makes it difficult to give a simple explanation of the 
SNR enhancement. The simplicity of the 1-dimensional Langevin equation allows us to 
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FIG. 5: Equilibrium distribution functions for the collective variable of a N = 10 units system 
with coupling parameter 9 = 0.5 and several D values obtained from the approximate Langevin 
equation Eq. ©. 

understand the reduction of the fluctuations in terms of the drastic differences between the 
drift coefficient D^(s) and D^(s,t) in driven systems. For driven systems, Eq. ([5]) also 
leads to a good approximation to the correlation function as seen in Fig. ([6]) where we depict 
the results for the incoherent part (upper panel) and the coherent part (lower panel) of the 
correlation function for S(t) as obtained from the simulation of Eq. ([1]) (solid lines) and 
the approximate Langevin equation, Eqs. (jSJ) and (J4]) (dashed lines) for N = 10, 6 = 0.5, 
D = 0.24 and a driving dichotomic force with amplitude A = 0.1 and fundamental frequency 
f2 = 0.01. is well approximated by (EI) while D^'(s, t) is fitted by fifth order polynomials 
for each half-period. 

In conclusion, we have analyzed a single collective variable characterizing a finite set 
of noisy bistable units with global coupling. We find several regions in parameter space 
separated by transition lines. The shape of P eq (s) switches from monomodal to multimodal 
as we move across the transition line. There are relevant qualitative differences with the 
results obtained in the infinite size limit, even though the lines separating different regions in 
parameter space look quite similar. We have also found approximate 1-dimensional Langevin 
equations for S(t). The change in the shape of P eq (s) implies the change on the drift 
coefficient as the parameter values are modified. In the presence of driving, the corresponding 
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FIG. 6: (Color online) The coherent (lower panel) and incoherent (upper panel) parts of the 
correlation function of S(t) in a system with N = 10, 9 = 0.5, D = 0.24 driven by a dichotomic 
force with A = 0.1 and f2 = 0.01. The results depicted in black correspond to the simulations 
of Eq. (pTJ) , while the red dashed lines are the results obtained with the 1-dimensional Langevin 
equation in ([5]). 

1-dimensional Langevin equation provides a reliable tool to understand the enhancement of 
the stochastic resonance effect in arrays relative to the one observed in a single bistable unit. 
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